#######################################################
# RUN MAIN IV MODEL WITH JEWISH INSURGENTS (DV: PROP) #
#######################################################
# Author: Kasia Nalewajko
# First created: 3 December 2021
# Replicated: 17 June 2024

rm(list = ls())

# LOAD PACKAGES -----------------------------------------------------------

if (!require("dplyr")) install.packages("dplyr")
if (!require("estimatr")) install.packages("estimatr")
if (!require("modelsummary")) install.packages("modelsummary")
if (!require("fixest")) install.packages("fixest")
if (!require("sf")) install.packages("sf")
if (!require("spdep")) install.packages("spdep")

# LOAD DATA ---------------------------------------------------------------

load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/main.Rda")

# RUN MODELS --------------------------------------------------------------

# model with limited controls

PROPIV_w_jew_ins_lim <- fixest::feols(log(((all_jewish_victims+1)/(allocated_jpop+1)+1)) ~ log(pop1936+1) + log(allocated_jpop+1) + synagogues_dummy + log(collabos_antijew1000+1) + log(DHI_milipol_sum1942+1) + area_sqkm  + longitude + longitudesq + latitude + latitudesq | zone5 + as.factor(ar_name)
                                      | log((FFIFFCall*1000/pop1936)+1) ~ log(mpf1000+1),
                                      cluster = ~c(id_bureau, resregion2),
                                      data = main)

# extract additional statistics
finfo <- fitstat( PROPIV_w_jew_ins_lim, "ivf")
PROPIV_w_jew_ins_lim_fstat <- round(finfo[["ivf1::log((FFIFFCall * 1000/pop1936) + 1)"]][["stat"]])
PROPIV_w_jew_ins_lim_wh <- round(summary( PROPIV_w_jew_ins_lim)$iv_wh[["p"]], digits = 4)

# run Moran test and extract the statistic
cantons_sf <- sf::st_read("./data/other/maps/france1940/CANTONS_1940/CANTONS_1940.shp")
obsRemoved <- abs(PROPIV_w_jew_ins_lim[["obs_selection"]][["obsRemoved"]])
cantons_sf <- cantons_sf[-obsRemoved,]
spatial_weights_matrix <- poly2nb(pl = cantons_sf, row.names = cantons_sf$deppct, queen = TRUE)
spatial_weights_matrix <- nb2listw(spatial_weights_matrix, zero.policy = T) 
PROPIV_w_jew_ins_lim_res <- residuals(PROPIV_w_jew_ins_lim)
PROPIV_w_jew_ins_lim_moran <- moran.test(x = PROPIV_w_jew_ins_lim_res, listw = spatial_weights_matrix, zero.policy = T)
PROPIV_w_jew_ins_lim_moran <- round( PROPIV_w_jew_ins_lim_moran$estimate[[1]], d = 3)

# fully controlled model

PROPIV_w_jew_ins_controls <- fixest::feols(log(((all_jewish_victims+1)/(allocated_jpop+1)+1)) ~ log(pop1936+1) + log(allocated_jpop+1) + synagogues_dummy + log(collabos_antijew1000+1) + log(DHI_milipol_sum1942+1)  + catholic_churches + log(FRANCISTE_36_1_pct+1) + log(ActionFrancaise_19_pct+1) + turnout_36_1_pct + log(nuance_droite_36_1_pct+1) + log(nuance_centredroit_36_1_pct+1) + log(nuance_centregauche_36_1_pct+1) + log(nuance_gauche_36_1_pct+1) + log(nuance_extremegauche_36_1_pct+1) + area_sqkm + longitude + longitudesq + latitude + latitudesq | zone5 + as.factor(ar_name) #  + log(nuance_extremedroite_36_1_pct+1)
                                           | log((FFIFFCall*1000/pop1936)+1) ~ log(mpf1000+1),
                                           cluster = ~c(id_bureau, resregion2), 
                                           data = main)

# extract additional statistics
finfo <- fitstat(PROPIV_w_jew_ins_controls, "ivf")
PROPIV_w_jew_ins_controls_fstat <- round(finfo[["ivf1::log((FFIFFCall * 1000/pop1936) + 1)"]][["stat"]])
PROPIV_w_jew_ins_controls_wh <- round(summary(PROPIV_w_jew_ins_controls)$iv_wh[["p"]], digits = 4)

# run Moran test and extract the statistic
cantons_sf <- sf::st_read("./data/other/maps/france1940/CANTONS_1940/CANTONS_1940.shp")
obsRemoved <- abs(PROPIV_w_jew_ins_controls[["obs_selection"]][["obsRemoved"]])
cantons_sf <- cantons_sf[-obsRemoved,]
spatial_weights_matrix <- poly2nb(pl = cantons_sf, row.names = cantons_sf$deppct, queen = T)
spatial_weights_matrix <- nb2listw(spatial_weights_matrix, zero.policy = T) 
mcontrols_county_res <- residuals(PROPIV_w_jew_ins_controls)
PROPIV_w_jew_ins_controls_moran <- moran.test(x = mcontrols_county_res, listw = spatial_weights_matrix, zero.policy = T)
PROPIV_w_jew_ins_controls_moran <- round(PROPIV_w_jew_ins_controls_moran$estimate[[1]], d = 3)

# PREVIEW -------

msummary(list(
  'First Stage' = PROPIV_w_jew_ins_lim$iv_first_stage[[1]],
  'Second Stage' = PROPIV_w_jew_ins_lim,
  'First Stage' = PROPIV_w_jew_ins_controls$iv_first_stage[[1]],
  'Second Stage' = PROPIV_w_jew_ins_controls
),
stars = c('*' = .1, '**' = .05, '***' = .01))

# EXPORT --------

# add_rows
add_stats <- tibble::tribble(
  ~x, ~y, ~q, ~w, ~l,
  'F stat. (1st stage)', PROPIV_w_jew_ins_lim_fstat, PROPIV_w_jew_ins_lim_fstat, PROPIV_w_jew_ins_controls_fstat, PROPIV_w_jew_ins_controls_fstat, 
  'Moran stat.', PROPIV_w_jew_ins_lim_moran, PROPIV_w_jew_ins_lim_moran, PROPIV_w_jew_ins_controls_moran, PROPIV_w_jew_ins_controls_moran,
  'Wu-Hausman p-value', PROPIV_w_jew_ins_lim_wh, PROPIV_w_jew_ins_lim_wh, PROPIV_w_jew_ins_controls_wh, PROPIV_w_jew_ins_controls_wh
)

attr(add_stats, 'position') <- c(48, 49, 50)

msummary(list(
  'First Stage' = PROPIV_w_jew_ins_lim$iv_first_stage[[1]],
  'Second Stage' = PROPIV_w_jew_ins_lim,
  'First Stage' = PROPIV_w_jew_ins_controls$iv_first_stage[[1]],
  'Second Stage' = PROPIV_w_jew_ins_controls
),
stars = c('*' = .1, '**' = .05, '***' = .01),
gof_omit = "AIC|BIC|RMSE|R2 Within|R2 Within Adj.",
coef_omit = "zone5",
add_rows = add_stats,
output = "latex",
coef_rename = c("log(mpf1000 + 1)" = "WWI military death rates",
                "fit_log((FFIFFCall * 1000/pop1936) + 1)" = "Insurgent presence (incl. Jewish insurgents)",
                "months_Verdun_Petain" = "Months Vérdun Pétain",
                "log(pop1936 + 1)" = "1936 population",
                "log(collabos_antijew1000 + 1)" = "Collaborators",
                "log(DHI_milipol_sum1942 + 1)" = "1942 state presence",
                "synagogues_dummy" = "Synagogues",
                "area_sqkm" = "Area size (km2)",
                "longitude" = "Longitude",
                "latitude" = "Latitude",
                "longitudesq" = "Longitude (sq)",
                "latitudesq" = "Latitude (sq)",
                "log(allocated_jpop + 1)" = "1941 Jewish population",
                "catholic_churches" = "Catholic churches",
                "log(FRANCISTE_36_1_pct + 1)" = "Franciste vote 1936",
                "log(ActionFrancaise_19_pct + 1)" = "Action Française vote 1919",
                "turnout_36_1_pct" = "Turnout 1936",
                "log(nuance_droite_36_1_pct + 1)" = "Right vote 1936",
                "log(nuance_centredroit_36_1_pct + 1)" = "Centre-right vote 1936",
                "log(nuance_centregauche_36_1_pct + 1)" = "Centre-left vote 1936",
                "log(nuance_gauche_36_1_pct + 1)" = "Left vote 1936",
                "log(nuance_extremegauche_36_1_pct + 1)" = "Extreme left vote 1936"
),
title = "Main results. Dependent variable: Proportion Holocaust victims (logged)"
)
